## Event Study Analysis (Results are estimated by lfe package version 2.8-7)
## These results are reported in Figure 4 and A2
rm(list=ls())
options(scipen=100,digits=10)
library(raster)
library(plyr)
library(dplyr)
library(lfe)
library(ggplot2)
library(stargazer)
load("data_final.RData")

########################
# Generate sample without post-shutdown period
data4=subset(data3, T<99)
########################
# Generate the short sample
data5=subset(data4, T>40)

##############################################################################################################################
## Event study for the short sample 
########################
## set reference week as the week before shutdown (redefine the group just in case...)
data5$Group="Control"
data5$Group[which(data5$Year==2018)]="Treated"
data5$Group[which(data5$weeks==9)]="Control"

########################
## Create interaction terms for estimating weekly gap coefficients
data5$interact=data5$Group
data5$interact[which(data5$Group=="Treated")]=paste(data5$Group[which(data5$Group=="Treated")], data5$weeks[which(data5$Group=="Treated")], sep = ":")

########################
## Event study regression for AOD
common_prod1=felm(AOD~  ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. + factor(weekday) + factor(Year) + factor(interact)|factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data5)
common=summary(common_prod1)
df=common_prod1$df
## Plot the results
alpha=as.data.frame(common$coefficients)
alpha$name=rownames(alpha)
interact=alpha[grepl(":",alpha$name),]
interact$low=interact$Estimate - qt(0.025,df)*interact$`Cluster s.e.`
interact$high=interact$Estimate + qt(0.025,df)*interact$`Cluster s.e.`
interact$Group="Insignificant"
interact$Group[which(interact$low<0&interact$high<0)]="Significant"
interact$Group[which(interact$low>0&interact$high>0)]="Significant"
interact$DATE=as.numeric(sub("\\:.*", "", gsub(".*?([0-9]).*?", "\\1", interact$name)))
ref=c(0,0,0,0,0,0,0,"Insignificant",9)
names(ref)=names(interact)
interact=rbind(interact, ref)
interact$DATE=as.numeric(interact$DATE)
interact$Estimate=as.numeric(interact$Estimate)
interact$low=as.numeric(interact$low)
interact$high=as.numeric(interact$high)
setwd("/Users/rjz5261/Dropbox/ShutDown/Final Submission Folder")
wd=getwd()
file.save.path <- file.path(wd,"Figures_Updated", "Event_AOD_short.jpg")
jpeg(file=file.save.path,
     width=1200,height=500)
ggplot(interact, aes(x = DATE, y = Estimate,color=Group),fill = "white") +
  geom_point(size = 5) +
  geom_errorbar(aes(ymax = high, ymin = low),width=0.3,linewidth=2)+
  scale_color_manual(values = c("grey55","black"))+
  geom_vline(xintercept=9, linetype="dashed", color='black', alpha=0.6)+
  annotate("text", x=9.7, y=0.1, label="EPA furlough", size=8, color='black', alpha=0.6)+
  geom_hline(yintercept=0) +
  scale_y_continuous(name="Weekly Differences between Treatment and Control Groups") +
  scale_x_continuous(breaks=c(6, 8, 10, 12, 14),
                     labels=c("Week of Dec. 2nd",  " Week of Dec. 16th",  "Week of Dec. 30th", "Week of Jan. 13th", "Week of Jan. 27th")) +
  theme(legend.position="bottom", panel.background=element_blank(),
        panel.border = element_rect(fill = NA, color = "black"),
        legend.text = element_text(size = 20),legend.key.size = unit(3,"line"),
        axis.text = element_text(size=15), axis.title = element_text(size=15)) +
  guides(color=guide_legend(""))

dev.off()

########################
## Event study regression for SO2
common_prod1=felm(SO2_MASS..tons.~  ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. + factor(weekday) + factor(Year) + factor(interact)|factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data5)
common=summary(common_prod1)
df=common_prod1$df
## Plot the results
alpha=as.data.frame(common$coefficients)
alpha$name=rownames(alpha)
interact=alpha[grepl(":",alpha$name),]
interact$low=interact$Estimate - qt(0.025,df)*interact$`Cluster s.e.`
interact$high=interact$Estimate + qt(0.025,df)*interact$`Cluster s.e.`
interact$Group="Insignificant"
interact$Group[which(interact$low<0&interact$high<0)]="Significant"
interact$Group[which(interact$low>0&interact$high>0)]="Significant"
interact$DATE=as.numeric(sub("\\:.*", "", gsub(".*?([0-9]).*?", "\\1", interact$name)))
ref=c(0,0,0,0,0,0,0,"Insignificant",9)
names(ref)=names(interact)
interact=rbind(interact, ref)
interact$DATE=as.numeric(interact$DATE)
interact$Estimate=as.numeric(interact$Estimate)
interact$low=as.numeric(interact$low)
interact$high=as.numeric(interact$high)
setwd("/Users/rjz5261/Dropbox/ShutDown/Final Submission Folder")
wd=getwd()
file.save.path <- file.path(wd,"Figures_Updated", "Event_SO2_short.jpg")
jpeg(file=file.save.path,
     width=1200,height=500)
ggplot(interact, aes(x = DATE, y = Estimate,color=Group),fill = "white") +
  geom_point(size = 5) +
  geom_errorbar(aes(ymax = high, ymin = low),width=0.3,linewidth=2)+
  scale_color_manual(values = c("grey55","black"))+
  geom_vline(xintercept=9, linetype="dashed", color='black', alpha=0.6)+
  annotate("text", x=9.7, y=2, label="EPA furlough", size=8, color='black', alpha=0.6)+
  geom_hline(yintercept=0) +
  scale_y_continuous(name="Weekly Differences between Treatment and Control Groups") +
  scale_x_continuous(breaks=c(6, 8, 10, 12, 14),
                     labels=c("Week of Dec. 2nd",  " Week of Dec. 16th",  "Week of Dec. 30th", "Week of Jan. 13th", "Week of Jan. 27th")) +
  theme(legend.position="bottom", panel.background=element_blank(),
        panel.border = element_rect(fill = NA, color = "black"),
        legend.text = element_text(size = 20),legend.key.size = unit(3,"line"),
        axis.text = element_text(size=15), axis.title = element_text(size=15)) +
  guides(color=guide_legend(""))

dev.off()

########################
## Event study regression for NOx
common_prod1=felm(NOX_MASS..tons.~  ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. + factor(weekday) + factor(Year) + factor(interact)|factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data5)
common=summary(common_prod1)
df=common_prod1$df
## Plot the results
alpha=as.data.frame(common$coefficients)
alpha$name=rownames(alpha)
interact=alpha[grepl(":",alpha$name),]
interact$low=interact$Estimate - qt(0.025,df)*interact$`Cluster s.e.`
interact$high=interact$Estimate + qt(0.025,df)*interact$`Cluster s.e.`
interact$Group="Insignificant"
interact$Group[which(interact$low<0&interact$high<0)]="Significant"
interact$Group[which(interact$low>0&interact$high>0)]="Significant"
interact$DATE=as.numeric(sub("\\:.*", "", gsub(".*?([0-9]).*?", "\\1", interact$name)))
ref=c(0,0,0,0,0,0,0,"Insignificant",9)
names(ref)=names(interact)
interact=rbind(interact, ref)
interact$DATE=as.numeric(interact$DATE)
interact$Estimate=as.numeric(interact$Estimate)
interact$low=as.numeric(interact$low)
interact$high=as.numeric(interact$high)
setwd("/Users/rjz5261/Dropbox/ShutDown/Final Submission Folder")
wd=getwd()
file.save.path <- file.path(wd,"Figures_Updated", "Event_NOx_short.jpg")
jpeg(file=file.save.path,
     width=1200,height=500)
ggplot(interact, aes(x = DATE, y = Estimate,color=Group),fill = "white") +
  geom_point(size = 5) +
  geom_errorbar(aes(ymax = high, ymin = low),width=0.3,linewidth=2)+
  scale_color_manual(values = c("grey55","black"))+
  geom_vline(xintercept=9, linetype="dashed", color='black', alpha=0.6)+
  annotate("text", x=9.7, y=1, label="EPA furlough", size=8, color='black', alpha=0.6)+
  geom_hline(yintercept=0) +
  scale_y_continuous(name="Weekly Differences between Treatment and Control Groups") +
  scale_x_continuous(breaks=c(6, 8, 10, 12, 14),
                     labels=c("Week of Dec. 2nd",  " Week of Dec. 16th",  "Week of Dec. 30th", "Week of Jan. 13th", "Week of Jan. 27th")) +
  theme(legend.position="bottom", panel.background=element_blank(),
        panel.border = element_rect(fill = NA, color = "black"),
        legend.text = element_text(size = 20),legend.key.size = unit(3,"line"),
        axis.text = element_text(size=15), axis.title = element_text(size=15)) +
  guides(color=guide_legend(""))

dev.off()



##############################################################################################################################
## Event study for the long sample 
########################
## set reference week as the week before shutdown (redefine the group just in case...)
data5=data4
data5$Group="Control"
data5$Group[which(data5$Year==2018)]="Treated"
data5$Group[which(data5$weeks==9)]="Control"
########################
## Create interaction terms for estimating weekly gap coefficients
data5$interact=data5$Group
data5$interact[which(data5$Group=="Treated")]=paste(data5$Group[which(data5$Group=="Treated")], data5$weeks[which(data5$Group=="Treated")], sep = ":")

########################
## Event study regression for AOD
common_prod1=felm(AOD~  ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. + factor(weekday) + factor(Year) + factor(interact)|factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data5)
common=summary(common_prod1)
df=common_prod1$df
## Plot the results
alpha=as.data.frame(common$coefficients)
alpha$name=rownames(alpha)
interact=alpha[grepl(":",alpha$name),]
interact$low=interact$Estimate - qt(0.025,df)*interact$`Cluster s.e.`
interact$high=interact$Estimate + qt(0.025,df)*interact$`Cluster s.e.`
interact$Group="Insignificant"
interact$Group[which(interact$low<0&interact$high<0)]="Significant"
interact$Group[which(interact$low>0&interact$high>0)]="Significant"
interact$DATE=as.numeric(sub("\\:.*", "", gsub(".*?([0-9]).*?", "\\1", interact$name)))
ref=c(0,0,0,0,0,0,0,"Insignificant",9)
names(ref)=names(interact)
interact=rbind(interact, ref)
interact$DATE=as.numeric(interact$DATE)
interact$Estimate=as.numeric(interact$Estimate)
interact$low=as.numeric(interact$low)
interact$high=as.numeric(interact$high)
setwd("/Users/rjz5261/Dropbox/ShutDown/Final Submission Folder")
wd=getwd()
file.save.path <- file.path(wd,"Figures_Updated", "Event_AOD_long.jpg")
jpeg(file=file.save.path,
     width=1200,height=500)
ggplot(interact, aes(x = DATE, y = Estimate,color=Group),fill = "white") +
  geom_point(size = 5) +
  geom_errorbar(aes(ymax = high, ymin = low),width=0.3,linewidth=2)+
  scale_color_manual(values = c("grey55","black"))+
  geom_vline(xintercept=9, linetype="dashed", color='black', alpha=0.6)+
  annotate("text", x=10, y=0.1, label="EPA furlough", size=8, color='black', alpha=0.6)+
  geom_hline(yintercept=0) +
  scale_y_continuous(name="Weekly Differences between Treatment and Control Groups") +
  scale_x_continuous(breaks=c(2, 4, 6, 8, 10, 12, 14),
                     labels=c("Week of Nov. 4th", "Week of Nov. 18th", "Week of Dec. 2nd",  " Week of Dec. 16th",  "Week of Dec. 30th", "Week of Jan. 13th", "Week of Jan. 27th")) +
  theme(legend.position="bottom", panel.background=element_blank(),
        panel.border = element_rect(fill = NA, color = "black"),
        legend.text = element_text(size = 20),legend.key.size = unit(3,"line"),
        axis.text = element_text(size=15), axis.title = element_text(size=15)) +
  guides(color=guide_legend(""))

dev.off()

########################
## Event study regression for SO2
common_prod1=felm(SO2_MASS..tons.~  ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. + factor(weekday) + factor(Year) + factor(interact)|factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data5)
common=summary(common_prod1)
df=common_prod1$df
## Plot the results
alpha=as.data.frame(common$coefficients)
alpha$name=rownames(alpha)
interact=alpha[grepl(":",alpha$name),]
interact$low=interact$Estimate - qt(0.025,df)*interact$`Cluster s.e.`
interact$high=interact$Estimate + qt(0.025,df)*interact$`Cluster s.e.`
interact$Group="Insignificant"
interact$Group[which(interact$low<0&interact$high<0)]="Significant"
interact$Group[which(interact$low>0&interact$high>0)]="Significant"
interact$DATE=as.numeric(sub("\\:.*", "", gsub(".*?([0-9]).*?", "\\1", interact$name)))
ref=c(0,0,0,0,0,0,0,"Insignificant",9)
names(ref)=names(interact)
interact=rbind(interact, ref)
interact$DATE=as.numeric(interact$DATE)
interact$Estimate=as.numeric(interact$Estimate)
interact$low=as.numeric(interact$low)
interact$high=as.numeric(interact$high)
setwd("/Users/rjz5261/Dropbox/ShutDown/Final Submission Folder")
wd=getwd()
file.save.path <- file.path(wd,"Figures_Updated", "Event_SO2_long.jpg")
jpeg(file=file.save.path,
     width=1200,height=500)
ggplot(interact, aes(x = DATE, y = Estimate,color=Group),fill = "white") +
  geom_point(size = 5) +
  geom_errorbar(aes(ymax = high, ymin = low),width=0.3,linewidth=2)+
  scale_color_manual(values = c("grey55","black"))+
  geom_vline(xintercept=9, linetype="dashed", color='black', alpha=0.6)+
  annotate("text", x=10, y=3.5, label="EPA furlough", size=8, color='black', alpha=0.6)+
  geom_hline(yintercept=0) +
  scale_y_continuous(name="Weekly Differences between Treatment and Control Groups") +
  scale_x_continuous(breaks=c(2, 4, 6, 8, 10, 12, 14),
                     labels=c("Week of Nov. 4th", "Week of Nov. 18th", "Week of Dec. 2nd",  " Week of Dec. 16th",  "Week of Dec. 30th", "Week of Jan. 13th", "Week of Jan. 27th")) +
  theme(legend.position="bottom", panel.background=element_blank(),
        panel.border = element_rect(fill = NA, color = "black"),
        legend.text = element_text(size = 20),legend.key.size = unit(3,"line"),
        axis.text = element_text(size=15), axis.title = element_text(size=15)) +
  guides(color=guide_legend(""))

dev.off()

########################
## Event study regression for NOx
common_prod1=felm(NOX_MASS..tons.~  ppt + tmean + tdmean + Wind + k_MWH+ SLOAD..1000.lbs. + factor(weekday) + factor(Year) + factor(interact)|factor(weeks)*factor(Plant_Code) + factor(T)|0|Plant_Code,
                  data = data5)
common=summary(common_prod1)
df=common_prod1$df
## Plot the results
alpha=as.data.frame(common$coefficients)
alpha$name=rownames(alpha)
interact=alpha[grepl(":",alpha$name),]
interact$low=interact$Estimate - qt(0.025,df)*interact$`Cluster s.e.`
interact$high=interact$Estimate + qt(0.025,df)*interact$`Cluster s.e.`
interact$Group="Insignificant"
interact$Group[which(interact$low<0&interact$high<0)]="Significant"
interact$Group[which(interact$low>0&interact$high>0)]="Significant"
interact$DATE=as.numeric(sub("\\:.*", "", gsub(".*?([0-9]).*?", "\\1", interact$name)))
ref=c(0,0,0,0,0,0,0,"Insignificant",9)
names(ref)=names(interact)
interact=rbind(interact, ref)
interact$DATE=as.numeric(interact$DATE)
interact$Estimate=as.numeric(interact$Estimate)
interact$low=as.numeric(interact$low)
interact$high=as.numeric(interact$high)
setwd("/Users/rjz5261/Dropbox/ShutDown/Final Submission Folder")
wd=getwd()
file.save.path <- file.path(wd,"Figures_Updated", "Event_NOx_long.jpg")
jpeg(file=file.save.path,
     width=1200,height=500)
ggplot(interact, aes(x = DATE, y = Estimate,color=Group),fill = "white") +
  geom_point(size = 5) +
  geom_errorbar(aes(ymax = high, ymin = low),width=0.3,linewidth=2)+
  scale_color_manual(values = c("grey55","black"))+
  geom_vline(xintercept=9, linetype="dashed", color='black', alpha=0.6)+
  annotate("text", x=10, y=1.5, label="EPA furlough", size=8, color='black', alpha=0.6)+
  geom_hline(yintercept=0) +
  scale_y_continuous(name="Weekly Differences between Treatment and Control Groups") +
  scale_x_continuous(breaks=c(2, 4, 6, 8, 10, 12, 14),
                     labels=c("Week of Nov. 4th", "Week of Nov. 18th", "Week of Dec. 2nd",  " Week of Dec. 16th",  "Week of Dec. 30th", "Week of Jan. 13th", "Week of Jan. 27th")) +
  theme(legend.position="bottom", panel.background=element_blank(),
        panel.border = element_rect(fill = NA, color = "black"),
        legend.text = element_text(size = 20),legend.key.size = unit(3,"line"),
        axis.text = element_text(size=15), axis.title = element_text(size=15)) +
  guides(color=guide_legend(""))

dev.off()